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Abstract: Turon 's methodology’ for determining optimal analysis parameters for the simulation of 
progressive delamination is reviewed. Recommended procedures for determining analysis 
parameters for efficient delamination growth predictions using the Abaqus/Standard cohesive 
element and relatively coarse meshes are provided for single and mixed-mode loading. The 
Abaqus cohesive element, COH3D8, and a user-defined cohesive element are used to develop 
finite element models of the double cantilever beam specimen, the end-notched flexure specimen, 
and the mixed-mode bending specimen to simulate progressive delamination growth in Mode I, 
Mode II, and mixed-mode fracture, respectively. The predicted responses are compared with their 
analytical solutions. The results show that for single-mode fracture, the predicted responses 
obtained with the Abaqus cohesive element correlate well with the analytical solutions. For 
mixed-mode fracture, it was found that the response predicted using COH3D8 elements depends 
on the damage evolution criterion that is used. The energy-based criterion overpredicts the peak 
loads and load-deflection response. The results predicted using a tabulated form of the BK 
criterion correlate well with the analytical solution and with the results predicted with the user- 
written element. 

Keywords: Damage, Composites, Progressive Delamination, Virtual Crack Closure Technique 
(VCCT), Cohesive Elements, Fracture Mechanics. 


1. Introduction 


Delaminations in laminated composite structures usually initiate from discontinuities such as 
matrix cracks and free edges or from embedded defects due to the manufacturing processes. 
Delaminations can reduce the stiffness and strength of a composite structure. In other cases, 
delamination can provide stress relief and delay the final failure of the composite structure. 
Therefore, it is important to analyze the progressive growth of delamination in order to predict the 
performance of a composite structure and to develop reliable and safe designs. 
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Abaqus/Standard provides the Virtual Crack Closure Technique (VCCT) and a cohesive element 
for the simulation of progressive delamination growth. The formulation of the cohesive element is 
based on the Cohesive Zone Model (CZM) approach for modeling complex fracture mechanisms 
at the crack tip (Dugdale, 1960 and Barenblatt, 1962). The cohesive zone approach models an 
extended cohesive zone, or process zone, ahead of the crack-tip using traction-separation laws that 
relate the opening displacements in the process zone to the resisting tractions. The traction- 
separation laws are defined in each fracture mode by an initial elastic stiffness, the peak traction, 
or interfacial strength, and the area under the traction-separation law that is equal to the critical 
energy release rate. 

Fracture simulations using cohesive elements require significant experience for determining the 
mesh requirements and the selection of properties that characterize the traction-separation 
response. For typical graphite-epoxy materials, the length of the cohesive zone is less than 1 mm. 
For an accurate representation of the distribution of tractions ahead of the crack-tip and an 
accurate simulation of delamination propagation, at least 3 to 5 cohesive elements are required in 
the cohesive zone. These mesh requirements may not be practical for the simulation of 
delamination propagation in structural analyses. 

Corigliano, 2003, used two-dimensional models of the Double Cantilever Beam specimen (DCB) 
and identified the pathological dependence of the solution on mesh refinement and attributed the 
errors in the predictions to the difficulties in calculating the interlaminar stresses accurately. Barua 
et al, 2007, conducted parametric studies of the effect of interfacial strength and mesh refinement 
on the predicted load-deflection response. Barua proposed establishing mesh requirements by 
performing mesh calibration studies using a simple DCB specimen and then applying those 
calibration results to more complex components made of the same material. 

Turon, 2007, proposed a procedure for relaxing the requirement for extremely fine meshes. The 
procedure consists of increasing artificially the length of the cohesive zone by reducing the 
interfacial strength. A longer cohesive zone allows for the use of larger elements while 
maintaining sufficient accuracy in the calculation of the energy release rate. 

Turon’ s procedure for using relatively coarse meshes has been demonstrated with a user-defined 
cohesive element (UEL) developed by Turon, 2006. In the present paper, Turon’s procedure is 
reviewed and is applied with the Abaqus cohesive element (COH3D8) to provide basic guidelines 
regarding optimal choices for mesh density and cohesive zone parameters for efficient single- and 
mixed-mode fracture simulations. Predictions for the load-displacement response for single and 
mixed-mode fracture specimens obtained using the Abaqus/Standard cohesive element are 
compared to the analytical solutions. 


2. Selection of Constitutive Parameters and Mesh Requirements 


The constitutive response of the cohesive elements used in this paper is defined by bi-linear 
traction-separation law as shown in Figure 1, where T represents the interfacial strength. A' is the 
critical separation, A fal1 is the separation at failure, and the area under the curve, G c , is the critical 
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strain energy release rate. Under mixed-mode fracture, the properties required to define the 
bilinear traction-separation law are the three critical fracture energies G IC , G IIC , G mc - the penalty 
stiffnesses Ki.Kj.Kt, and the interfacial strengths T, Si, and ,S\ In the present work, the shear 
strengths in the two orthogonal directions, S ) and S 2 are assumed to be the same and are referred to 
as S. 

The initial response of the cohesive element is assumed to be linear until a damage initiation 
criterion is met. The penalty stiffness, K h of the bi-linear traction-separation law is defined as 


where A'- is the critical separation for damage initiation. The value of the penalty stiffness must be 
high enough to prevent interpenetration of the crack faces and to prevent artificial compliance 
from being introduced into the model by the cohesive elements. However, an overly high value 
can lead to numerical problems. Several guidelines have been proposed for obtaining the penalty 
stiffness of a cohesive element. Turon, et al., 2007, proposed an equation for the penalty stiffness 
for Mode I that ensures that the compliance of the bulk material is much larger than the initial 
compliance of the cohesive element: 


K 3 >a ~ (2) 

where a is a parameter much larger than 1, E 3 is the transverse Young’s modulus of the material, 
and t is the thickness of an adjacent sub-laminate. Turon recommends using a = 50 . This value is 
used for the penalty stiffness in all of the analyses presented herein. In the present paper, the 
penalty stiffnesses for all modes are assumed to be the same: K_=K 2 =K 3 =K. 



Figure 1. Bilinear traction-separation response of the cohesive element. 

To predict the propagation of delamination accurately, it is necessary to have a sufficiently fine 
finite element discretization of the cohesive zone to ensure a good representation of the 
interlaminar stress fields. The required maximum size of the cohesive elements can be determined 
by evaluating the size of the cohesive zone and the number N e of elements that are needed in the 
zone. 
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The length of the cohesive zone, l cz , is the distance from the crack tip to the point when the 
maximum cohesive traction is attained. Figure 2 describes the length of the cohesive zone. It has 
been shown that, for infinite-sized geometries, the length of the cohesive zone, l cz , is a material 
property. For single Mode I and Mode II fracture, the cohesive lengths are approximately (Turon, 
2007) 


Lj = ME 2 and l czJI = ME 2 ^L 


(3) 


where E 2 is the transverse Young’s modulus and M is a parameter that depends on the cohesive 
zone theory used to derive the expression for the cohesive zone length. In this study, Hillerborg’s 
model in which M is equal to 1.0 is used (Hillerborg, 1976). Therefore, the length I e of the 
cohesive elements must satisfy the following conditions for Mode I and Mode II: 


l<- 


L 


N„ 


1 and l< lcz 11 


N„ 


(4) 



F,A 


F, A 


Figure 2. Length of the cohesive zone. 

Turon suggests that a minimum of three elements in the cohesive zone is required to accurately 
represent the fracture energy in the cohesive zone (Turon, 2007). However, the cohesive zone 
length for typical polymeric composites is less than 1 mm, and such small mesh requirements may 
render most structural analyses computationally intractable. 

When the length of the cohesive zone is negligible compared to the crack length, the propagation 
of damage is governed by linear elastic fracture mechanics (LEFM) principles that state that the 
energy release rate for creating a new fracture area must be equal to G c . Therefore, when LEFM is 
valid, the effect of the interfacial strengths can be neglected. It can be observed in Equation 3 that 
the length of the cohesive zone increases as the inverse of the square of the interfacial strength. 
This suggests that the length of the cohesive zone can be artificially lengthened to ensure that it 
spans enough elements of a given size. The interfacial normal and shear strengths that are required 
can be obtained using Equations 3 and 4, which give 
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T 


^2 ' G IC 


and S' 


^2 ' Gnc 


( 5 ) 


N-L 


N.-l. 


Therefore, the normal and shear strengths that ensure that enough elements span the cohesive zone 
are 


T =Min{T a j) and 5 =Min{s a ,S } (6) 

Abaqus provides four damage initiation criteria: maximum stress, maximum strain, quadratic 
stress, and quadratic strain criterion. The quadratic stress criterion, shown in Equation 7, is used 
throughout this paper. 



( 7 ) 


where the Macaulay bracket, < >, signifies that the compressive stress does not contribute to 
damage initiation. 

Once the damage initiation criterion of the cohesive element is satisfied, the stiffness of the 
element is degraded. Equation 8 represents the softening response of the cohesive element 

<Tj = (1 - d)K i A i , i = 1,2,3 (8) 


where d is the damage variable, which has the value d= 0 when the interface is undamaged, and 
the value d = 1 when the interface is fully fractured. 

In the analyses presented herein, we use the energy-based Benzeggagh and Kenane (BK) damage 
evolution criterion shown in Equation 9 (Benzeggagh and Kenane, 1996) 


G C =G IC + {G 1IC -G ll 


J Shear 


( 9 ) 


where // is the BK material parameter, G shear = G u + G m , and G r = G , + G n + G m . 

In the UEL used herein (Turon, 2006), the BK fracture criterion can be expressed in terms of 
separations: 


Gc=G ic +{G U c-G 


Ic ) 


p 2 

1 + 2 P 2 -ip 


where /? is the mode mixity ratio defined in terms of separation. 


(10) 
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p= 


A 


shear 


A 


shear 


+ <A 


normal 


) 


(ii) 


Abaqus provides two methods to implement the BK criterion. In the standard BK option, the 
fracture toughness G c given by Equation 9 is determined in terms of the accumulated energy 
release rates G shear and G r . This option is specified using the following command lines 

♦DAMAGE EVLOUTION, TYPE=ENERGY, 

MIXED MODE BEHAVIOR=BK, POWER= rj 

Alternatively, the BK criterion can also be specified with G c supplied in tabular form. Using 
Equation 10, a table of G c can be specified as a function of the mode mixity ratio. This option is 
specified using the following command lines 

♦DAMAGE EVLOUTION, TYPE=ENERGY, 

MIXED MODE BEHAVIOR=Tabular 

G c , <h 

where is defined as 


<k=- tan-’ (Ml-/*)) (12) 

K 

In the following simulations, both options are used and their corresponding results are compared 
to the results obtained with a user-defined element and with analytical solutions. 

3. Analyses of progressive delamination growth 


The Abaqus cohesive element COH3D8 and the user-defined cohesive element COHUEL 
developed by Turon, 2006 were used to develop finite element (FE) models of the double 
cantilever beam (DCB) specimen and the end-notched flexure (ENF) specimen to simulate 
progressive delamination growth under Mode I and Mode II fracture, respectively. In Sections 3.1 
and 3.2, the results of nonlinear analyses performed to investigate the effect of mesh density and 
interfacial strength on the predicted load-deflection response of the DCB and ENF specimens are 
presented. From the progressive delamination growth predictions of the FE models for single- 
mode fracture, the mesh density required for solution convergence and the corresponding normal 
and shear strengths for a given material are obtained. Progressive delamination growth predictions 
of single-mode fracture obtained with the Abaqus cohesive element and COH UEL correlate 
closely with each other. Therefore, only the progressive delamination growth predictions obtained 
with the Abaqus cohesive element are presented with the analytical solutions. 

FE models of the mixed-mode bending (MMB) specimen were developed using the mesh density 
and the interfacial strengths that resulted in converged results for the single-mode analyses. In 
Section 3.3, the predicted progressive delamination results for the MMB specimen with the 
Abaqus cohesive element and the COH UEL are compared with analytical solutions and with 
Abaqus VCCT predictions. 
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FE models of each specimen are developed using the material properties of AS4/3501-6 that are 
provided in Table 1. All specimens are unidirectional and the fibers are aligned with the direction 
of fracture propagation. 


Table 1. Mechanical and interface material properties of AS4/3501-6. 


E, 

e 2 

e 3 

V 12 

Gl2 

Gl3 

G 23 

148 GPa 

10.50 GPa 

10.50 GPa 

0.27 

5.61 GPa 

5.61 GPa 

3.17 GPa 


Gic 

Giic 

Giiic 

T 

S 

n 

0.08 kJ/nf 

0.55 kJ/rrf 

0.55 kJ/rrf 

53.78 MPa 

86.88 MPa 

1.8 


3.1 Double Cantilever Beam Specimen 

FE models of a DCB specimen are developed to investigate the influence of mesh refinement and 
the normal interfacial strength on the predicted responses for Mode I simulation. The length of the 
specimen, L, is 101.6 mm long, and the width is 7.62 mm. The specimen is composed of two sub- 
laminates of thickness t=l .524 mm. The initial crack length, a 0 , is 29.21 mm. Each sub-laminate 
has a stacking sequence of [0 3 ]. The configuration of the specimen is shown in Figure 3. 

The DCB specimen is modeled with two layers of 4-node shell elements S4 that are connected 
together with cohesive elements. The reference surfaces of the sub-laminates are moved from the 
mid-surface of the sub-laminates to the contact surface between the two sub-laminates using the 
shell property option OFFSET 

Three levels of mesh refinement were considered and the corresponding models are labeled DCB- 
1, -2, and -3. All three models have a uniform mesh in a 10 mm-long damage propagation zone 
ahead of the crack tip, which is labeled c 0 in Figure 3. The element sizes in zone c Q of models 
DCB-1, DCB-2, and DCB-3 are 1.0 mm, 0.5 mm, and 0.32 mm, respectively. 


1 Crack tip 

z 



1 



^ a 0 c 0 w \ 

r L * 


Figure 3. Configuration of double cantilever beam (DCB) specimen. 

The cohesive zone length for Mode I fracture calculated from Equation 3 is l cz /= 0.29 mm. The 
penalty stiffness determined from Equation 2 with a = 50 is K = 344 GPa/mm. The nominal normal 
strength provided in Table 1 and the adjusted normal strength obtained using Equation 5 are used 
to study the influence of mesh size and normal strength on the prediction of fracture propagation 
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in Mode I. Three different values for the number of elements within the cohesive zone, N e , are 
substituted into Equation 5 to calculate the adjusted normal strength. The number of elements 
within the cohesive zone and the adjusted normal strength corresponding to the different mesh 
refinements, DCB-1, -2, and -3, are shown in Table 2. 

Table 2. Adjusted normal strengths corresponding to the three mesh refinements. 


Number of Elements in 
Cohesive Zone 

DCB-1 

le Model ~ 1 mm 

DCB-2 

le model ~ 0.5 mm 

DCB-3 

le model ~ 0.32mm 

N e -1 

29.6 MPa 

41.9 MPa 

53.8 MPa 

N e =3 

17.1 MPa 

24.2 MPa 

30.2 MPa 

N e =5 

13.2 MPa 

18.7 MPa 

23.4 MPa 


The predicted load-deflection responses obtained using model DCB-1 are compared to the 
analytical solution (Reeder, 2004) in Figure 4. The model with the nominal normal strength and an 
adjusted normal strength with N e =1 significantly overpredicts the strength of the DCB specimen 
and the response does not follow the softening response of the analytical solution. The responses 
predicted using an adjusted strength with N e =3 and 5 compare well with the analytical solution. 
However, these models exhibit over-softening before the peak load and they display spurious 
oscillations during delamination propagation. Therefore, model DCB-1 is too coarse to predict 
accurately the propagation of delamination, even when the strengths are lowered using Equation 5. 

The predicted load-deflection responses obtained using model DCB-2 are shown in Figure 5. The 
response predicted with the nominal normal strength and the response predicted with an adjusted 
normal strength with N e = 1 significantly overpredict the peak load. Both responses exhibit minor 
oscillations during fracture propagation. The responses and the peak loads obtained with N e = 3 
and N e — 5 correlate well with the analytical solution. 

The load-deflection responses obtained using the most refined model, DCB-3, are shown in Figure 
6. For this level of mesh refinement, the adjusted strength with N e =1 is approximately equal to the 
nominal strength. The solution obtained with nominal strength overpredicts the peak load and the 
softening response of the analytical solution by approximately 7%. The peak loads and the load- 
deflection responses predicted using normal strengths adjusted with N e = 3 and N e = 5 correlate 
well with the analytical solution. 

A comparison of the results obtained with models DCB-1, DCB-2, and DCB-3 indicates that, 
when used in conjunction with an adjusted strength with N e = 3, model DCB-2 is the coarsest 
model that provides accurate results. 
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Figure 4. Effect of the normal interfacial strength on the predicted load-deflection 

response for FE model DCB-1. 



Figure 5. Effect of the normal interfacial strength on the predicted load-deflection 

response for FE model DCB-2. 
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Figure 6. Effect of the normal interfacial strength on the predicted load-deflection 

response for FE model DCB-3. 


3.2 End Notched Flexure (ENF) Specimen 


FE models of the ENF specimen were developed to investigate the influence of mesh refinement 
and the interfacial strength on the predicted response for Mode II fracture simulation. The ENF 
specimen has the same dimensions and stacking sequence as the DCB specimen analyzed in 
Section 3.1. For the ENF specimen, the force, F, is applied at the middle of the specimen. The 
ENF specimen configuration and the boundary conditions are shown in Figure 7. 


FE models of the ENF specimens with two levels of mesh refinement, indicated as model ENF-1 
and ENF-2, are developed to investigate the sensitivity of the predicted responses to the mesh 
refinement. These models have a uniform mesh in a damage propagation zone of length c D equal to 
21.6 mm ahead of the crack tip. The span-wise element length in zone c 0 of Models ENF-1 and 
ENF-2 is 1.0 mm and 0.5 mm, respectively. The cohesive zone length for Mode II fracture 
calculated from Equation 3 is l cz n =0.77 mm. 


Crack tip 




w =0 


) 

X 


1 


0.5L 


ao 


Co 


2t 


u, w =0 


Figure 7. Configuration and boundary conditions of the End Notched Flexure 

specimen. 
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A comparison of the predicted load-deflection responses of the FE models and the analytical 
solution (Mi, 1998) is shown in Figure 8. The results of the FE models for the ENF specimen are 
consistent with the analytical solutions and indicate that both models are sufficiently refined and 
that the difference in element size between the two models has a minimal influence on the 
predicted responses. In addition, when using an element size of 1 mm or smaller, no adjustment to 
the nominal shear strength is necessary. These results indicate that the simulation of progressive 
delamination of a composite structure subjected to Mode II 1 does not require a mesh as fine as a 
mesh required for Mode I fracture. 



Figure 8. Effect of mesh refinement on the predicted load-deflection response of 

ENF FE models. 


3.3 Mixed-Mode Bending (MMB) Specimen 

The Abaqus cohesive element and the COH UEL are used to develop FE models of the MMB 
specimen. The specimen is 100 mm-long and 25.5 mm-wide with two 2.3 mm-thick sub- 
laminates. The initial crack length, a 0 , is 27 mm. The experimental set-up of Reeder, et al., 1990, 
shown in Figure 9 is modeled. The level of mode mixity in the MMB is applied by varying the 
length of the applied load lever arm, d, of the test apparatus. 

In the FE simulation of the MMB specimen, the lever arm length, d, is equal to 52.8 mm and the 
corresponding mixed-mode ratio, G U /G T , is equal to 0.2. The penalty stiffness determined from 
Equation 2 with a=50 and t=2.3 mm is equal to A=228GPa/mm. 
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The FE model of the MMB specimen has a uniform mesh in a damage propagation zone ahead of 
the crack tip of length c 0 equal to 23.0 mm. The element size in the propagation zone was selected 
to be l e = 0.5 mm, and the number of elements in the cohesive zone is chosen as N e = 3. These 
values were found to provide converged results in the analyses of the DCB and ENF specimens 
described in the previous sections. Models of the MMB were created with the Abaqus cohesive 
elements and with user-defined elements. Two damage evolution options were evaluated with the 
Abaqus cohesive elements: the standard BK and the tabular options described in Section 2. To 
improve the convergence rate, a viscoelastic regularization factor of 10" 4 was used in the analyses 
for both the Abaqus cohesive element and the user-defined element. 

Using Equation 5, the reduction factor of the normal strength, f IT , is 45 %. The ENF analyses 
indicate that the nominal shear strength of the material does not need to be reduced. Flowever, to 
investigate the influence of the shear strength on the mixed-mode fracture predictions, three 
different shear strengths are used: the nominal shear strength, S , the shear strength adjusted by 
using Equation 5, S , and the shear strength reduced by the same reduction factor as the normal 
strength, S. The adjusted normal and shear strengths are reported in Table 3. 


Table 3. The shear strengths of the cohesive element. 



AS4/3501-6 

T : Adjusted Normal Strength by Equation 5 

24.2 MPa 

S : Nominal Shear Strength 

86.9 MPa 

S : Shear Strength Adjusted by Equation 5 

63.7 MPa 

S : Shear Strength Adjusted by Mode 1 reduction factor 

39.2 MPa 



L 


Figure 9. Configuration, experimental set-up, and boundary conditions of the MMB 

specimen. 

The predicted load-deflection responses obtained using Abaqus cohesive elements, the VCCT 
tool, and the analytical solution are shown in Figure 10. The VCCT response and the analytical 
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solution correlate well with each other. The responses predicted with the standard BK option and 
the Abaqus COH3D8 elements overpredict the force required for crack propagation by 
approximately 35%. However, when the tabulated form of the BK option is used, the predicted 
results correlate well with the analytical solution and with the VCCT prediction. The predicted 
solutions are also sensitive to the normal-to-shear strengths used in the analyses. For the ratios 
used, the three load-deflection curves obtained with the tabulated BK option differ by 
approximately 6%. 

The results predicted with the standard BK option and with the tabulated BK option should be 
identical to each other if the mode-mixity does not change during damage evolution. The fact that 
the results shown in Figure 10 do not correlate well with each other indicates that the mode-mixity 
changes during damage evolution. This unexpected result will require further investigation of the 
evolution of damage in the MMB test. 

The predicted load-deflection responses obtained using the user-defined cohesive element (Turon, 
2006), COHUEL, are shown in Figure 11. When the adjusted normal and nominal shear 
strengths, T, S, are used, the FE model overpredicts the peak load and the load-deflection response 
compared to the analytical solution and the VCCT prediction. The results also indicate that 
reducing the shear strength influences the predicted load-deflection response of the FE model. 
When the adjusted normal and shear strengths, f , S, obtained by Equation 5 are used for the 
interfacial strengths, the predicted load-deflection responses obtained with the COH UEL 
correlate well with the analytical solution and with Abaqus VCCT predictions. Although the 
nominal shear strength can be used as shear strength of COH UEL for single Mode II fracture, the 
predicted response of FE models with COH UEL indicates that both normal and shear strengths 
needed to be adjusted using Turon’s guideline to obtain accurate fracture simulations for the 
mixed-mode fracture. 



Figure 10. Load-deflection response of MMB specimen with Abaqus COH3D8 

cohesive element. 
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Figure 11. Load-deflection response of MMB specimen with user-defined COH_UEL 

cohesive element. 


4. Concluding Remarks 


Turon’s methodology for selecting analysis parameters for the simulation of delamination 
propagation using relatively coarse meshes is reviewed and used to determine analysis parameters 
for use with the Abaqus/Standard cohesive element. The Abaqus cohesive element, COH3D8, and 
a user-defined cohesive element are used to model a DCB specimen, an ENF specimen, and a 
MMB specimen to simulate progressive delamination in Mode I, Mode II, and mixed-mode 
fracture, respectively. For single-mode fracture, the predicted responses obtained with the Abaqus 
cohesive elements correlate well with the analytical solutions when the interfacial normal and 
shear strengths are determined using Turon’s methodology. For mixed-mode fracture, the 
predicted response obtained with Abaqus cohesive element depends on the damage evolution 
criterion selected. The results using the standard BK overpredict the analytical results by as much 
as 35%, but when the BK criterion is supplied in tabulated form, the results correlate well with the 
VCCT results and with the analytical solution. 

The selection of cohesive element parameters and mesh requirements for fracture propagation in 
structural subcomponents should be established by performing single-mode fracture simulations in 
simple specimens of the same material. In the present paper, the optimal mesh density and analysis 
parameters for the simulation of fracture of a mixed-mode bending test were determined from 
parametric simulations of DCB and ENF specimens. This procedure shoidd also be applied for the 
simulation of fracture in more complex structures. 
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